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ABSTRACT 

We present general relativistic solutions for self-similar spherical perturbations in an 
expanding cosmological background of cold pressure-less gas. We focus on solutions 
having shock discontinuities propagating in the surrounding cold gas. The pressure, 
p, and energy-density, /i, in the shock-heated matter are assumed to obey p — wfi, 
where w is a positive constant. Consistent solutions are found for shocks propagating 
from the symmetry center of a region of a positive density excess over the background. 
In these solutions, shocks exist outside the radius marking the event horizon of the 
black hole which would be present in a shock-less collapse. For large jumps in the 
energy-density at the shock, a black hole is avoided altogether and the solutions are 
regular at the center. The shock-heated gas does not contain any sonic points, provided 
the motion of the cold gas ahead of the shock deviates significantly from the Hubble 
flow. For shocks propagating in the uniform background, sonic points always appear 
for small jumps in the energy-density. We also discuss self-similar solutions without 
shocks in fluids with w < —1/3. 
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1 INTRODUCTION 



Detailed studies of generic self-gravitating systems are only 
possible with the aid of complex numerical methods. Part of 
the complication arises from the long range nature of grav- 
itational interactions. Thus, even equilibrium states, when 
they exist, do not in general obey standard statistical me- 
chanics. This is correct for both Newtonian gravity as well 
as general relativity (GR), but the latter adds its pecu- 
liar technical and conceptual challenges. Analytic or semi- 
analytic treatment of systems possessing special symme- 
tries are, therefore, of great interest. Gravity is the dom- 
inant force dictating the growth of structure in the uni- 
verse. Thus, self-gravitating perturbations in an expanding 
cosmological has traditionally gained a great deal of atten- 
tion (cf. Peebles 1990). Here we present a class of solutions 
which describe the evolution of spherical perturbations in 
an expanding Einstein-de Sitter cosmological background. 
We focus on the evolution of self-similar perturbations in 
GR. Self-similarity occurs when dimensionless physical vari- 
ables can be expressed as functions of some combination 
of the time and space coordinates. Self-similar systems of- 
fer great simplification since they are described by ordinary 
differential equations. Newtonian treatment of cosmological 
self-similar perturbations have been presented in great detail 



previously in the literature (e.g. Fillmore & Goldreich 1984; 
Bertschinger 1985; Chuzhoy & Nusser 2000). Nevertheless, 
it is instructive to review the general features of these solu- 
tions. In cosmology, we assume that the density fluctuations 
are extremely small at time, ti, close to the the Big Bang sin- 
gularity. An isolated small density perturbation having the 
form S/j, ~ r~ s (r is the distance from the symmetry cen- 
ter) at ti, develops a self-similar behavior at sufficiently late 
times, t ^> ti* . Consider first a positive density perturbation 
in a universe made solely of dust (i.e. collision-less matter). 
As time goes-by a shell of matter begins to expand at a rate 
lower than the background until it reaches a maximum dis- 
tance from the center. The shell then falls inward and ends 
up oscillating around the center together with shells which 
have reached it at earlier times. According to this picture, 
the radius, n, of the shell which is at maximum expansion 
at time t plays an important role. For r > n, matter is still 
expanding, while at r < ri shells of matter are oscillating 
in the central potential well generated by their self-gravity. 



* An additional assumption which is often made is that the initial 
expansion rate of matter at r is equal to that of the unperturbed 
background. This assumption is, however, unnecessary as long as 
the deviations from the background expansion rate are not great. 
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The self-similar variable can then be taken to be the ratio, 
r/ri(t). For a universe filled with a collisional fluid (gas), 
self-similarity develops provided the initial pressure is neg- 
ligible and that the gas is able to shock and evolve adia- 
batically. In this case, an inner region of hot shocked gas 
develops. The shock bounding this region propagates away 
from the center at a rate dri(t)/ dt. 

Self-similarity in GR is obtained for a more restricted 
class of initial perturbations than in Newtonian gravity. GR 
admits self-similar solutions only for specific values of the 
power law index s, while Newtonian dynamics does not im- 
pose such restrictions. Further, the pressure, p, must be 
related linearly to the energy-density, jj,, of matter (Cahil 
& Taub 1971). Carr & Hawking (1974) and Carr & Yahil 
(hereafter CY90) have studied self-similar perturbations in 
an expanding universe. Their treatment, however, did not 
include shocks which might form as a result of collapse or 
energy injection into the system. Here we complement the 
work of CY90 by exploring solutions which involve shock 
discontinuities. We assume that the isotropic cosmic back- 
ground is made of cold pressure-less matter that is able to 
shock. We will present solutions in which shocks propagate 
from a central point into the uniform cosmic background or 
into the outer regions of a perturbation with positive excess 
of energy-density in the cold gas. Shock-heated gas will be 
assumed to obey the equation of state p = wfi (w > 0). 
We also present solutions with no shocks for w < —1/3, a 
case which has not been discussed at length in the literature 
before. 

The paper is organized as follows. In §2 we summarize 
GR equations for self-similar perturbations in comoving co- 
ordinates. For a detailed discussion of these equations we 
refer the reader to CY90 and Carr & Coley (2000) (here- 
after CC02). In §3 we discuss the formation of shocks and 
the jump conditions relating the physical variables on both 
sides of the shock discontinuity. In §4 we present full nu- 
merical solutions. We conclude with a general discussion in 
§5. 



2 THE SELF-SIMILAR GR EQUATIONS IN 
SPHERICAL SYMMETRY 

We work with comoving coordinates as they are particularly 
convenient for the study of cosmological perturbations. With 
only minor changes of notation our description of the equa- 
tions closely follows CY90 and CC02. Denoting the time, 
the radius and the two spherical angular coordinates by t, 
r, 6 and <j>, the space-time metric is expressed in terms of 
the interval ds as (e.g. Landau & Lifshitz 1975), 



, 2 2k j, 2 

ds = e dt 



e 2X dr 2 - R 2 dQ 2 



(1) 



where v, A and R are, in general, functions of t and r, and 
dQ. 2 — d8 2 +sin 2 9 d(f> 2 ■ In comoving coordinates the energy- 
momentum tensor has the form T^ v = (fj, + p)u^u v — pg^ v 
where u M = dt/ ds — (e~ v , 0, 0, 0). For an equation of state 
p — w/j, where w is a constant, self-similar solutions in which 



dimensionless variables are functions of z = r/t alone can 
be obtained. We use the following dimensionless variables^, 

S(z) = | , £(*) = ( M r 2 r 1/(1+lu) and M(z) =^,(2) 
where 



m(r,t) = ^ 



(3) 



Using T$. a = 0, gives d t X + 2d t R/R = -2d t n/(p + n) and 
d r v = —2d r p/(p + fi). For p — w/j,, these equations yield 



and e A = 7^ 1 S 2 , 



(4) 



Substituting these relations in the remaining Einstein equa- 
tions and using d/dt = —(z/t)d/ dz and d/dr = (z/r) d/ dz 
yield 



s + s+[s + (i + w )s] (Y^f-f I 



( 2WJ 2 \ c 4 , 2 S e 2-2w (2-2iu)/(1+to) 



(5) 



1 + w 



-(V 2 -w) 7 5 4 | = (l + ro )^ 1 - t ", 



l + (l + w)- 



(6) 
(7) 



M = S 2 C (1+W) 

M = \ + _^ r 2»^(i-»)/d+») _ ±*fC 2 S° ( 1 + I ] («) 



and 



V = {(3 1 )- 1 ^ 1 - w S- 2 z {1 - w)/(1+w) , 



(9) 



where an over-dot denotes a derivative with respect to y = 
In z. 

A particular solution of these equations describes the 
Einstein-e Sitter universe, i.e. the fiat critical density 
Friedmann-Robertson- Walker (FRW) uniform density back- 
ground. This solution is 

£ = z -2/(l+w) ^ g = z -2/(3+3«0 



with 



V2 



_ 3 + 3w 
P ~ V3{l + w) ' 7 ~ l + 3w ' 
The space-time interval is 
ds 2 = p 2 dt 2 - 7 - 2 z- 4/(3+3w) dr 2 

_ r 2(l+3»)/(3+3») t 4/(3+3») dQ 2 

This is brought to a standard form 
ds 2 = dP - i i/(3+3w) ( df 2 + f 2 dfi 2 ) 
using the transformation 
t = fit and r = [3 



2/(3+3m) r (l+3u>)/(3+3u;) 



(10) 
(11) 

(12) 
(13) 
(14) 



t The variable £ defined here is related to x{z) in CY90 by £ 1 / v 
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2.1 Equations in perturbed variable 



3 SHOCKS 



In order to describe deviations from the isotropic FRW back- 
ground, we follow CY90 and introduce A(z) and B(z) ac- 
cording to"f 



j- _ -2/(l+w) A 

£ = z ' K ' e 



S = z 



-2/(3+3™) B 



(15) 



From the definition of £ in (2), we get fit 2 = e A . Thus, 
to linear order A can be identified as the energy-density 
constrast Sfi/fi. The equations (5) through (8), expressed in 
terms of A and B, are, 



B = - A + (1 + w)AB - 3B 2 



/ l + 3w \ 
\ 1 + w J 



B 



M 



- + {l + w)B 



e 2B-(l+w)A z 2(l+3w)/(3+3w) 



(16) 
(17) 
(18) 



We will consider solutions involving shock discontinuities 
propagating in a medium of cold pressure-less gas. We will 
assume that the heated matter behind the shock obeys the 
equation of state p — wfi. A shock will be specified by the 
ratio of the energy-densities on both of its sides. The re- 
maining hydrodynamical variables and the components of 
the metric on both sides of the shock are then related by 
means of the jump conditions. These are usually derived 
by working first with Schwarzschild coordinates where the 
metric tensor is continuous and then transforming back to 
comoving coordinates (Bardeen 1965; Cahil & Taub 1971). 
Quantities just ahead and behind the shock are denoted by 
the (+) and (— ) symbols, respectively. The general form of 
the jump conditions in comoving coordinates is (Cahil & 
Taub 1971), 



M 



_ 1 3(1 + w) 2 c 2B-2wA z 2(l+3w)/(3+3w) 



2 4 

\2 



B 



9(1 + w) 2 
2(1 + 3m) 2 

and 



B + 



l + 3w 



3(1 + 



V = 



l+3w \ (l+3«,)/(3+3«,) (1-w)A-2B 



3(1 + TO 

(20) 
(21) 



(19) 



2.2 Asymptotic solutions 

Linear perturbations describing small deviations from uni- 
formity are obtained for A ~ < 1 and B C 1. The 
linearized equations have two solutions. The first is (CY90) 



A 1 
Si 



1 + 3™ 



k z 



-2(1+3to)/(3+3uj) 



1+W 

B - fc z - 2 ( 1 + 3 ™)/( 3 + 3 ™) 



where the constants Boa and k are related by 

, _ 3(l + »)(e- 2S ~ -e 4S ~) 
k ~ 2 



(l + 3w)(5 + 3w) 



(22) 
(23) 

(24) 



Expressed in terms of the FRW coordinates f and t given in 
(14), we see that Ai ~ Sfi/n oc f-2p(i+3«0/(3+3«0 _ There _ 
fore, the energy-density contrast, <5/u//u, is f~ 2 , independent 
of w, and has a temporal behavior identical to that obtained 
for large scale perturbations in the standard linear perturba- 
tion theory in the comoving gauge (e.g. Peebles 1993). The 
second solution is 

(i-«0/(i+«0 „„j D ^ _ ^2(2) 



^2(2) oc z 



and B2(z) 



(25) 



According to (14), the limit r — » 00 corresponds to 2 — > cm 
and 2 — > 0, respectively for w > —1/3 and — 1 < w < —1/3. 
Therefore, density perturbations approaching zero as f — > 00 
are described by Ai for w > —1/3 and A2 for — 1 < w < 
— 1/3. The mode Ai grows with time, while A2 decays. 

i Note that the A is equal to w times the function used by CY90. 



v+ 



V- 



V- = 



(W-fl- ~ W + (1 + )()I- + W+fJ,+ ) 

(H- - n+){p,+ +w-n-) 
(w-n- - w + fi + )( f i+ + w-n-) 

(fx- - fl+){fJ,- +W+fl+) 
(/£- + W + fl+)(fl+ + W+fJ,+ ) 
(fl+ + W-fl-)(p- + W-fl-) 

(ju+ + w-fi-)(n+ + w+n+) 
(fi- +w+fi+)(n- +w-fi-) 

c + r+ 

1 + CT+ 



where 



. v dR/dt 
dR/dr 



and C ■■ 



V+ - V- 
1 - V-V+ 



The condition (26) leads to the continuity of the self- 
variable S, i.e. S+ — S-. Using the relations 



(26) 
(27) 

(28) 

(29) 

(30) 

(31) 

(32) 
similar 



Rt = rS t = 
R r = S + rS r = S + r 
we find that 



dS dz ■ 
r— — = — zS , 



d2 dt 

d^dz 
dz dr 



S + S . 



r = 



-V(*)S 
S + S 



(33) 



Therefore, the condition (31) relates between S = 2 dS/ dz 
on both sides of the shock. 



4 NUMERICAL SOLUTIONS 

We present here numerical solutions of the equations of mo- 
tion for cases with and without shocks. The boundary condi- 
tions for all solutions obey the asymptotic form giving small 
density perturbations in the limit of large FRW coordinate 
f, corresponding to 2 > 1 mi > -1/3 and to 2 <JC 1 for 
w < -1/3. 
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Figure 1. Numerical solutions of self-similar positive density per- 
turbations in a FRW universe for w = (solid lines), -0.2 (dotted), 
-0.3 (dashed) and -0.6 (dash-dotted) . Top-left: the excess energy- 
density in units of the background value. Top-right: the function 
M = S/R. Bottom-left: the fluid "velocity" AR/ dt = -z 2 dS/dz. 
Bottom-right: the velocity V of surfaces of constant z relative to 
the fluid. For w = —0.6, the lines showing M, dR/ dt and V are 
scaled down by factors of 10 4 , 500 and 25, respectively. 



4.1 Solutions without shocks 

The case w > has been considered extensively by CY90. 
Therefore, we only present solutions for — 1 < w < 0. An im- 
portant difference between solutions with w > and w < 
is the following. For w > 0, a sonic point at which V = w" 1 ^ 2 
may exist. At this point the term (1 — w/V 2 ) in equation (17) 
vanishes, leaving A undetermined. Certain physical assump- 
tions must then be made in order to complete the solution 
across the sonic point (e.g. Ori & Piran 1990; CY90). For 
w < 0, sonic points do not exist and, barring singularities, 
A is always determined by (17). 

Figure 1 shows results of numerical integration of a pos- 
itive density perturbation for w — (solid lines), w = —0.2 
(dotted), w — —0.3 (dashed), and w = —0.6 (dash-dotted). 
The four panels show, from top-left in clockwise order, the 
contrast in the energy-density (i.e. fi/jl — 1 where p, is the 
background energy density), the function M — m/R, the 
"velocity" dR/ dt = -z 2 dS/dz, and the velocity V of 
surfaces of constant z relative to the fluid. Solutions with 
w > —1/3 diverge at a finite value of z. This singular be- 
havior of the solutions indicates the presence of black holes 
with apparent horizon at the point with M = 1/2 and an 
event horizon at the point with V = 1 (Carr & Hawking 
1974). Comoving coordinates should allow a description of 
the physical variables down to the singularity at the cen- 
ter. For the solutions with w > —1/3, the finite value of 
z at which the solutions diverge represents the shell hav- 
ing a comoving coordinate r and arriving at the center at 
time t — r/z (e.g. Landau & Lifshitz 1975). The lines of 
t 2 n - 1 follow the asymptotic form z (1 - w)/(1+w) at z < 1 



for w < -1/3 and 2 2(i+3»)/(3+3») at 2 > f for w > -1/3. 
For w < —1/3, the limit z — * is obtained by taking the 
limit f — > oo or/and t cx^ oo. Therefore, the self-similar 
solution for w = —0.6 describes a mode which decays with 
time and approaches zero at large FRW coordinate f, in 
contrast to the solutions shown for w > —1/3 which grow 
with time. In the Newtonian limit the quantity dR/ dt de- 
scribes the velocity of a shell relative to a central observer. 
For w > —1/3, at sufficiently large z, dR/ dt is positive. As 
we decrease z, dR/ dt changes signs, indicating the collapse 
of matter towards the center. 

4.2 Solutions with shocks 

We discuss now solutions describing self-similar perturba- 
tions in the presence of shocks. The cosmic background is 
assumed to be made of cold pressure-less gas which is able 
to shock. We focus first on shocks propagating from the cen- 
ter of a perturbation with positive energy-density excess in 
the cold gas. At the end of this section we present solutions 
for shocks propagating into the uniform density background. 
The equations of motion (16) and (17) with w = are inte- 
grated inward from zi 1 with initial conditions dictated 
by the linear solution. The integration gives tabulated val- 
ues for the variables A(z) and B(z) from zi until z = z s 
(z s < zi), where a shock is assumed to exist. Given the 
variables A+ and B+ just a head of the shock at z a we ob- 
tain the corresponding physical variables v + , A+, S+ and 
£+. Then, assuming a value for the ratio g = fj,-/fj, + of the 
energy-densities behind and ahead of the shock, we use the 
jump conditions (26) through (31) to obtain A_, S- 
and £_ behind the shock. These variables are then used as 
initial conditions in the numerical integration of the equa- 
tions (5) and (6), from z s down to smaller values. This is 
done for several values g and z s . The condition (27) implies 
that shocks cannot exist at points with V > 1. This is re- 
lated to that fact that V = 1 marks an event horizon inner 
to which matter must eventually hit the central singularity 
(Carr & Hawking 1974). Results of the numerical integra- 
tions are shown in figure (2). The top-left panel plots the 
energy-density (in units of the baclground value), t 2 /i. The 
dashed line represents the solution obtained for the shock- 
less collapse of cold gas (w — 0). The solid lines show t 2 fi 
versus z when shocks are assumed to exist at lnz a = —3.46 
and In z s — —4.05. For these z s , the shock-less solution has 
dR/ dt < (the bottom panel to the left). For each value 
of z s , lines representing solutions obtained for four values 
of g = n-/fi+ are shown: g — 1.1, 2, 3, and 5. At both z s , 
the energy- density reaches a constant value as z — > for 
the two larger values of g and diverges at a finite z for the 
smaller values. Given g, the jump conditions determine the 
pressure and energy-density ahead of the shock and, there- 
fore, the ratio w = W- = p_//n_. The top-right panel plots 
lines of w versus g for three values of y a — In z a as indi- 
cates in the figure. The behavior of these lines can be easily 
understood by a straightforward algebraic manipulation of 
the jump condition (27) with w + = 0. For g close to unity, 
w increases linearly with g — 1. In the limit 3 > 1 we get 
w w V+ 2 /(l — V+ 2 )/g so that p_ = wfx- which is indepen- 
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Figure 2. Results of numerical integration of the equations in 
the presence of shocks propagating from the center of a positive 
density excess. The dashed lines in the top-left, bottom-left and 
bottom-right panels correspond to the shock-less collapse of cold 
gas (u! = 0) . The solid lines in these panels represent the physical 
variables in the shock-heated region. Results are shown for shocks 

present at z s = —3.46 and z s = —4.05; for g = /u. = 1-2, 

2, 3, and 5 at each z s . For these values of z s the shock-less so- 
lution has dR/ dt < 0. Top-left: the energy-density in units of 
the background density versus the variable y = In 2. Top-left: 
the parameter w in the shocked material as a function of the as- 
sumed g, for shocks at three values of z s as indicated in the figure. 
Bottom-left: the "velocity" dR/ dt = —z 2 dS/ dz. Bottom-right: 
the velocity V of constant z surfaces relative to the fluid, where 
the dashed line is the actual velocity of the collapsing cold gas 
(w = 0), while the solid lines show V/w 1 / 2 . 

dent of g. Curves of dR/ dt are shown in the bottom-left 
panel. The numeric labels on each line indicate the corre- 
sponding value of g. It is interesting that for the regular 
solutions, dR/ dt converges to a positive constant as z — > 0. 
This constant becomes smaller with decreasing z s and we 
conjecture that it will approach zero as z s gets very close 
to the event horizon marked by the point with V = 1 in 
the dashed curve (bottom panel to the left). The physical 
interpretation of this is that outflows/winds from a central 
region are required for the formation of shocks. The inten- 
sity of these outflows determine the location of the shock, z s . 
The bottom-right panel plots the velocity V versus z. The 
actual velocity is shown for the cold collapse (dashed) while 
the ratio V/w 1 ^ 2 is plotted for the other lines. For the regular 
solutions, V/w 1 / 2 is less then unity throughout the shocked 
regions. The other solutions have V/w 1 ^ 2 > 1 over the range 
from z a down to the finite z marking the central singularity. 
Thus, in both cases no sonic points exist. Finally, the asymp- 
totic behavior of the numerical solutions near the center is 
consistent with the analysis of CC02. For the regular solu- 
tions, we have V — > as z — > 0. The analysis of CC02 of this 
case gives S oc z^ 1 and £ oc z~ 2 ^ 1+w ^ which according to 
(2) implies t 2 /i — const and dR/ dt = — z 2 dS/ dz = const. 
Our solutions match this behavior extremely well. 
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Figure 3. Comparison between solutions for shocks propagat- 
ing in the uniform cosmic background (bottom panels) and in a 
positive density excess over the background (top). In all panels, 
dashed lines represent solutions without shocks, while the solid 
ones correspond to shocks assumed to exist at two different posi- 
tions. The solid lines at each shock position correspond to several 
values of g = /j._//u_|_ from g = 1.1, to 6 in equal steps of 0.1. Left 
panels show lines of energy-density versus y = In 2. Right panels: 
show V/w 1 / 2 for the solid lines and the actual the velocity V for 
the dashed line. 

It is interesting to compare between the solutions for 
shocks propagating in a positive density perturbation and 
in the uniform background. These solutions are represened, 
respectively, in the top and bottom rows of in figure (??). 
The positive density solutions are taken from the output of 
the integration used in figure (2) but for shocks assumed 
to exist at two points outside the region with dR/ dt < 
(compare with figure (2). This means that the shock exist 
at a point where the cold fluid is expanding. Physically this 
is possible if a central energy source is present. As is seen in 
the top panel to the right, solid curves with intermediate g 
hit sonic points (V = to 1 / 2 ). At these points we stop the in- 
tegration. Larger values of g yield regular solutions without 
sonic points. For smaller g, the velocity V changes direc- 
tion and tends to infinity, indicating the presence of a black 
hole. As z s is decreased, the solutions behave in a manner 
similar to what we have seen in figure (2). In the uniform 
background case, sonic points always exist for small g and 
disappear for sufficiently large values. 



5 DISCUSSION 

In the standard cosmological model, the universe on large 
scales is nearly uniform so that super-horizon perturbations 
are well described by the general relativistic theory of linear 
perturbations (e.g. Bardeen 1980). However, in view of re- 
sent observational data, it seems that GR might be needed 
for a full study of perturbations on small scales where strictly 
speaking the linearized equations are invalid. These data 
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suggest that the expansion of the universe is accelerating 
(e.g. Spergel et. al. 2003). It is thought that this acceler- 
ation is driven by a dark-energy component which, in its 
simplest incarnation, is described in terms of a fluid with 
pressure p related to the energy density, /i, by p — w[i with 
w w — I at the present cosmic epoch. The continuity equa- 
tion, p — —(I + w)fi implies that for w = — 1 the density is 
constant with time and, therefore, by Lorentz invariance, the 
fluid must be homogeneous. More generally, dark energy is 
modeled in terms of a scalar field which, under certain con- 
ditions, could be represented as a fluid with w varying with 
time. Such a fluid is not ruled out by observations. If this 
is the case, inhomogencities in the dark energy component 
could develop. Since dark energy is inherently a relativis- 
tic entity, GR therefore should be invoked for a thorough 
treatment of small scale cosmological perturbations. Seen 
in this perspective self-similar solutions might be relevant 
for the description of structure formation in the universe. 
However, the solutions are quite limited in their applicabil- 
ity in a realistic model for the growth of structure in the 
presence of dark energy. Putting the subdominant baryonic 
matter aside, the problem of structure formation requires 
tracing the evolution of the gravitationally coupled system 
of collisionless dark matter and dark energy. Our solutions 
do not take this coupling into account. It is unclear whether 
self-similar evolution of cosmological perturbations in a two 
fluid system is possible at all. Further, we assume a strict 
equation of state of the form p — wp. This leads to an imag- 
inary "sound" speed, (Sp/Sfi) 1 ^ 2 for w < 0, while a more 
realistic model for dark energy in terms of scalar fields leads 
in general to a non-imaginary sound speed. Verertheless, be- 
cause of the paucity of analytic solutions for time-dependent 
systems in GR, a treatment of systems with restricted sym- 
metry is a worthwhile effort. Analytic study of these sys- 
tems could also be useful in testing and calibrating numeri- 
cal techniques for solving the GR equations. 
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